
clear all

cd "$revdta"

set more off

scalar define exclusion_restriction = 1		/*1=Baseline, 2=Simulated welfare benefits , 3=Simulated UI benefits, 4=No excl. restr. */ /*Changing to 2 or 3 produces Appendix Table 7; changing to 4 produces App Table 8 col. (4)*/ 
scalar define SRC_only 				= 0		/*0=Baseline, 1=Drop SEO*/ /*Changing to 1 produces Appendix Table 5*/

******************************************************************************************************************************************
use "Full_Panel_2b.dta", clear

keep if sample_analytic==1 & educ <= 2

drop yrd* age2
*NOTE: to get original results, just restrict to years 2008 and prior in this data.
*need to get rid of yrd18-21 too
*yes this does it largely, at least for reg coefs

merge m:1 newid using ".\intermediary\sp_temp_types"
drop lw_res
scalar def lowlim		=0.85
scalar def highlim		=5

if("$type" != "type_prodwage") ren type_prodwage type_prod
tokenize "$typename"
*"post"
local type $type  
*type_prodpost
ren `type' type_prodwage


/*SPOUSE - Drop those with outlier wage growth, below 1/2 of mw, or with less than 3 obs, or outside relevant age range*/
sort sp_newid year
qby sp_newid:gen grly=(exp(sp_lw)/exp(sp_lw[_n-1]))-1
replace grly=. if grly==-1
gen suspect=0
qui su grly,d

replace suspect=grly<-lowlim|(grly>highlim & grly!=.)	
egen ss=sum(suspect),by(newid)				/*# of outlier wage growth records*/
egen sw=sum(exp(sp_lw)<0.5*minwage),by(newid)	/*#of wages below state-level minimum wage*/
egen numy=sum(newid!=.),by(newid)			/*# of obs*/

gen sp_todrop=sw>0|numy<3|ss>0				/*Drop those with outlier wage growth, below 1/2 of mw, or with less than 3 obs*/

sort newid year

replace sp_lw = . if sp_todrop==1
replace sp_lw = . if agew < 23 | agew > 62
if exclusion_restriction==2	{
	replace all_benefits=all_wf
	}
else if exclusion_restriction==3	{
	replace all_benefits=ui_pay
	}
else if exclusion_restriction==4	{
	replace all_benefits=0
	}
	
if SRC_only==1	{
	drop if seo==1
	}

	*PLOT EARNINGS YOU'D HAVE IF WORKING FOR REAL
	
gen sp_part    =sp_lw!=. if sp_todrop != 1 & agew >= 23 & agew <= 62
gen part    =lw!=.

gen sp_sev = sp_DS==2
gen sp_mod = sp_DS==1 | sp_DS==2
gen sp_dis = sp_sev==1|sp_mod==1

g all_sev=all_benefits*sev
g all_mod=all_benefits*mod
g all_welf=all_benefits*(year>=1996)						
g all_sev_welf=all_benefits*(year>=1996)*sev
g all_mod_welf=all_benefits*(year>=1996)*mod
g all_rec=all_benefits*(year>=2008)						
g all_sev_rec=all_benefits*(year>=2008)*sev
g all_mod_rec=all_benefits*(year>=2008)*mod


g sp_all_sev=all_benefits*sp_sev
g sp_all_mod=all_benefits*sp_mod
g sp_all_dis=all_benefits*sp_dis
g sp_all_welf = all_benefits*(year>=1996)
g sp_all_dis_welf=all_benefits*(year>=1996)*sp_dis
g sp_all_sev_welf=all_benefits*(year>=1996)*sp_sev
g sp_all_mod_welf=all_benefits*(year>=1996)*sp_mod
g sp_all_rec = all_benefits*(year>=2008)
g sp_all_dis_rec=all_benefits*(year>=2008)*sp_dis
g sp_all_sev_rec=all_benefits*(year>=2008)*sp_sev
g sp_all_mod_rec=all_benefits*(year>=2008)*sp_mod



****************************************DEFINE 0/1 ESTIMATION METHOD ***************
tab year,g(yrd)
g age2=age^2/100	
replace agew=. if agew==999 | agew==0
g sp_age = agew
g sp_age2 = agew^2/100

sort sp_newid year
foreach var of varlist sp_lw sp_sev sp_mod sp_dis  sp_age sp_age2 {
qby sp_newid: g diff`var'= `var'-`var'[_n-1]
}

sort newid year
qby newid: g gap=year-year[_n-1]
foreach var of varlist lw sev mod age age2 white married yrd1-yrd21 {
	qby newid: g diff`var'=`var'-`var'[_n-1]
	replace diff`var' = . if !(gap==1 | gap == 2)
	}

tab type_prodwage, gen(prodtype_)

g ageh = age*(mod==0)*(sev==0)
g sp_ageh = sp_age*(sp_dis==0)*(sev==0)
g agemod = age*mod
g sp_agemod = sp_age*sp_dis
g agesev= age*sev

g age2h = age2*(mod==0)*(sev==0)
g sp_age2h = sp_age2*(sp_dis==0)*(sev==0)
g age2mod = age2*mod
g sp_age2mod = sp_age2*sp_dis
g age2sev= age2*sev

g old = age >= 45

g agehold = age 
replace agehold = 0 if age >= 45 | sev == 1 | mod == 1
g age2hold = age2 
replace age2hold = 0 if age >= 45 | sev == 1 | mod == 1

g sp_agehold = sp_age
replace sp_agehold = 0 if sp_age >= 45 | sp_dis == 1
g sp_age2hold = sp_age2 
replace sp_age2hold = 0 if sp_age >= 45 | sp_dis == 1

g agemodold = age 
replace agemodold = 0 if age >= 45 | mod == 0
g age2modold = age2 
replace age2modold = 0 if age >= 45 | mod == 0 

g agesevold = age 
replace agesevold = 0 if age >= 45 | sev == 0
g age2sevold = age2 
replace age2sevold = 0 if age >= 45 | sev == 0 

g sp_agemodold = sp_age 
replace sp_agemodold = 0 if sp_age >= 45 | sp_dis == 0
g sp_age2modold = sp_age2 
replace sp_age2modold = 0 if sp_age >= 45 | sp_dis == 0 


g hold = age >= 45 & sev == 0 & mod == 0
g diffsev_old = 0
replace diffsev_old = diffsev if age >= 45
g sev_old = 0
replace sev_old = sev if age >= 45
g diffmod_old = 0
replace diffmod_old = diffmod if age >= 45
g mod_old = 0
replace mod_old = mod if age >= 45

g sp_old = sp_age >= 45

g sp_hold = sp_age >= 45 & sp_dis==0
g diffsp_mod_old = 0
replace diffsp_mod_old = diffsp_dis if age >= 45
g sp_mod_old = 0
replace sp_mod_old = sp_dis if age >= 45


g vold = age >= 50
g diffsev_vold = 0
replace diffsev_vold = diffsev if age >= 50
g diffmod_vold = 0
replace diffmod_vold = diffmod if age >= 50
g sev_vold = 0
replace sev_vold = sev if age >= 50
g mod_vold = 0
replace mod_vold = mod if age >= 50

g vvold = age >= 55
g diffsev_vvold = 0
replace diffsev_vvold = diffsev if age >= 55
g diffmod_vvold = 0
replace diffmod_vvold = diffmod if age >= 55
g sev_vvold = 0
replace sev_vvold = sev if age >= 55
g mod_vvold = 0
replace mod_vvold = mod if age >= 55

******************************************************************************************************************************************
**********PARTICIPATION DECISION 1/0- Table 2 column 1
*selection correction: benefits interacted with disability and pre-post 1996 welfare reform
forvalues m=0(1)1{
dprobit part sev mod age age2 white all_benefits all_sev all_mod all_welf all_sev_welf all_mod_welf yrd1-yrd21 if married==`m' , vce(cluster newid) 
if exclusion_restriction!=4	{
	testparm all*
}
predict xb0_m`m',xb
g lamda0_m`m'=normalden(xb0_m`m')/normal(xb0_m`m') if married==`m'
}
*overall head:
dprobit part sev mod age age2 white married all_benefits all_sev all_mod all_welf all_sev_welf all_mod_welf yrd1-yrd21 , vce(cluster newid) 
predict xb0,xb
g lamda0=normalden(xb0)/normal(xb0)

*spousal:
dprobit sp_part sp_dis sp_age sp_age2 white all_benefits sp_all_dis sp_all_welf sp_all_dis_welf yrd1-yrd21 if married==1 , vce(cluster newid) 
predict sp_xb0,xb
g sp_lamda0=normalden(sp_xb0)/normal(sp_xb0) if married==1

***SELECTION CORRECTION WITH POST-2008 DUMMY
*overall head:
dprobit part sev mod age age2 white married all_benefits all_sev all_mod all_welf all_sev_welf all_mod_welf all_rec all_sev_rec all_mod_rec yrd1-yrd21 , vce(cluster newid) 
predict xb0_rec,xb
g lamda0_rec=normalden(xb0_rec)/normal(xb0_rec)

*spousal:
dprobit sp_part sp_dis sp_age sp_age2 white all_benefits sp_all_dis sp_all_welf sp_all_dis_welf sp_all_rec sp_all_dis_rec yrd1-yrd21 if married==1 , vce(cluster newid) 
predict sp_xb0_rec,xb
g sp_lamda0_rec=normalden(sp_xb0_rec)/normal(sp_xb0_rec) if married==1



gen bothpart=part*sp_part

preserve
collapse (mean) part, by(mod sev age)
graph twoway connected part age, by(sev mod)
restore
preserve
collapse (mean) sp_part, by(sp_dis sp_age)
graph twoway connected sp_part sp_age if sp_age >= 23 & sp_age <= 62, by(sp_dis)
restore

**********WITH GROUP FE AND AGE JUMP INSTRUMENT XXX****
*overall head:
probit part prodtype_* sev mod ageh age2h agehold age2hold agemod age2mod agemodold age2modold agesev age2sev agesevold age2sevold white married yrd1-yrd21  all_benefits all_sev all_mod all_welf all_sev_welf all_mod_welf hold mod_old sev_old, vce(cluster newid) nocons
eststo headwork
predict xb0_fe,xb
g lamda0_fe=normalden(xb0_fe)/normal(xb0_fe)
testparm all* hold mod_old sev_old

quietly probit part prodtype_* sev mod ageh age2h agehold age2hold agemod age2mod agemodold age2modold agesev age2sev agesevold age2sevold white married yrd1-yrd21  all_benefits all_sev all_mod all_welf all_sev_welf all_mod_welf hold mod_old sev_old if sev == 1, vce(cluster newid) nocons
margins, dydx(sev_old) at(sev=1)

**********WITH GROUP FE****
*by spouse
forvalues m=0(1)1{
dprobit part prodtype_* sev mod age age2 white all_benefits all_sev all_mod all_welf all_sev_welf all_mod_welf yrd1-yrd21 if married==`m' , vce(cluster newid) 
if exclusion_restriction!=4	{
	testparm all*
}
predict xb0_fem`m',xb
g lamda0_fem`m'=normalden(xb0_fem`m')/normal(xb0_fem`m') if married==`m'
}

*spousal:
ren old oldhead
ren hold holdhead
ren mod_old mod_oldhead
ren mod modhead
ren age agehead
ren age2 age2head
ren ageh agehhead
ren age2h age2hhead
ren agehold ageholdhead
ren age2hold age2holdhead
ren agemod agemodhead
ren age2mod age2modhead
ren agemodold agemodoldhead
ren age2modold age2modoldhead
ren all_sev all_sevhead 
ren all_mod all_modhead 
ren all_sev_welf all_sev_welfhead
ren all_mod_welf all_mod_welfhead
ren all_welf all_welfhead

ren sp_old old
ren sp_hold hold
ren sp_mod_old mod_old
ren sp_dis mod
ren sp_age age
ren sp_age2 age2
ren sp_ageh ageh
ren sp_age2h age2h
ren sp_agehold agehold
ren sp_age2hold age2hold
ren sp_agemod agemod
ren sp_age2mod age2mod
ren sp_agemodold agemodold
ren sp_age2modold age2modold
ren sp_all_dis all_mod
ren sp_all_welf all_welf
ren sp_all_dis_welf all_mod_welf


**********SPOUSAL WITH GROUP FE AND AGE JUMP INSTRUMENT XXX****
*spouse:
probit sp_part prodtype_*  mod ageh age2h agehold age2hold agemod age2mod agemodold age2modold white yrd1-yrd21 all_benefits all_mod all_welf all_mod_welf hold mod_old if married==1 , vce(cluster newid) nocons
eststo spousework

********************************************************

estimates restore spousework
predict sp_xb0_fe,xb
g sp_lamda0_fe=normalden(sp_xb0_fe)/normal(sp_xb0_fe) if married==1
testparm all*

ren old sp_old
ren hold sp_hold 
ren mod_old sp_mod_old 
ren mod sp_dis 
ren age sp_age 
ren age2 sp_age2 
ren ageh sp_ageh 
ren age2h sp_age2h 
ren agehold sp_agehold 
ren age2hold sp_age2hold 
ren agemod sp_agemod 
ren age2mod sp_age2mod 
ren agemodold sp_agemodold 
ren age2modold sp_age2modold 
ren all_mod sp_all_dis 
ren all_welf sp_all_welf 
ren all_mod_welf sp_all_dis_welf 
ren all_welfhead all_welf 

ren oldhead old
ren holdhead hold 
ren mod_oldhead mod_old 
ren modhead mod 
ren agehead age 
ren age2head age2 
ren agehhead ageh 
ren age2hhead age2h 
ren ageholdhead agehold 
ren age2holdhead age2hold 
ren agemodhead agemod 
ren age2modhead age2mod 
ren agemodoldhead agemodold 
ren age2modoldhead age2modold 
ren all_sevhead  all_sev 
ren all_modhead  all_mod 
ren all_sev_welfhead all_sev_welf 
ren all_mod_welfhead all_mod_welf 



*joint: 
dprobit bothpart prodtype_* sev mod ageh age2h agehold age2hold agemod age2mod agemodold age2modold agesev age2sev agesevold age2sevold sp_dis sp_ageh sp_age2h sp_agehold sp_age2hold sp_agemod sp_age2mod sp_agemodold sp_age2modold white all_benefits  all_sev all_mod all_welf all_sev_welf all_mod_welf hold mod_old sev_old sp_all_dis sp_all_welf sp_all_dis_welf  sp_hold sp_mod_old  yrd1-yrd21 if married==1 , vce(cluster newid) 
predict both_xb0_fe,xb
g both_lamda0_fe=normalden(both_xb0_fe)/normal(both_xb0_fe) if married==1
**************************
cd "$out"


esttab headwork spousework using "workparams`1'_jump.tex", tex fragment keep(all*  hold mod_old sev_old mod sev age*  married prodtype_*)  coeflabels(sev "Sev. Dis." mod "Mod. Dis." ageh "Age " age2h " $ Age^2/100$" married "Married" prodtype_1 "Low Type" prodtype_2 "Medium Type" prodtype_3 "High Type") replace mtitles("Head" "Partner") nonum se eqlabels(, none)
estout headwork using "workparams`1'_jump.txt", keep(all*  hold mod_old sev_old mod sev age*  married prodtype_*) replace nonum eqlabels(, none)  cells("b se p")  mlabels(,none) 
estout spousework using "spworkparams`1'_jump.txt",  keep(all*  hold mod_old mod age*  prodtype_*) replace nonum eqlabels(, none)  cells("b se p")  mlabels(,none) 

sort newid year
/*
forvalues m=0(1)1{
qby newid: g lag1lamda0_m`m'=lamda0_m`m'[_n-1]
qby newid: g lag1xb0_m`m'=xb0_m`m'[_n-1]
}
*/
qby newid: g lag1lamda0=lamda0[_n-1]
qby newid: g lag1xb0=xb0[_n-1]

qby newid: g lag1sp_lamda0=sp_lamda0[_n-1]
qby newid: g lag1sp_xb0=sp_xb0[_n-1]

qby newid: g lag1lamda0_rec=lamda0_rec[_n-1]
qby newid: g lag1xb0_rec=xb0_rec[_n-1]

qby newid: g lag1sp_lamda0_rec=sp_lamda0_rec[_n-1]
qby newid: g lag1sp_xb0_rec=sp_xb0_rec[_n-1]

qby newid: g lag1lamda0_fe=lamda0_fe[_n-1]
qby newid: g lag1xb0_fe=xb0_fe[_n-1]

forvalues m=0(1)1{
qby newid: g lag1lamda0_fem`m'=lamda0_fem`m'[_n-1]
qby newid: g lag1xb0_fem`m'=xb0_fem`m'[_n-1]
}

qby newid: g both_lag1lamda0_fe=both_lamda0_fe[_n-1]
qby newid: g both_lag1xb0_fe=both_xb0_fe[_n-1]


qby newid: g lag1sp_lamda0_fe=sp_lamda0_fe[_n-1]
qby newid: g lag1sp_xb0_fe=sp_xb0_fe[_n-1]

********FD, WITH SELECTION CORRECTION - Table 2 column 2
*interacts the probit fit and its lags with time gap 
*(to treat pre-post change in PSID differently)
*why lags of probit fit? stickiness of work decision?
g dyreff=0
g yreff=0
g dyreff_nosel=0
g yreff_nosel=0
g dyreff_rec=0
g yreff_rec = 0
g dyreff_fe=0
g yreff_fe = 0
g yreff_m1=0
g dyreff_m1=0
g yreff_m0=0
g dyreff_m0=0
g sp_yreff=0
g sp_g1=.
g sp_dyreff=0
g sp_reslev=.
g sp_yreff_nosel=0
g sp_g1_nosel=.
g sp_dyreff_nosel=0
g sp_reslev_nosel=.
g sp_yreff_rec=0
g sp_g1_rec = .
g sp_dyreff_rec=0
g sp_reslev_rec = .
g sp_yreff_fe=0
g sp_g1_fe = .
g sp_dyreff_fe=0
g sp_reslev_fe = .


*Overall Wage equation for the household head, no selection:
xi: reg difflw diffsev diffmod diffage diffage2   diffmarried diffwhite diffyrd* ,vce(cluster newid) nocons
eststo headresults_noselection

replace dyreff_nosel=. if !e(sample)
replace yreff_nosel=. if !e(sample)
forvalues j=1(1)21	{
	replace dyreff_nosel=dyreff_nosel+_b[diffyrd`j']*diffyrd`j' if e(sample)
					}
gen g1_nosel=difflw-_b[diffsev]*diffsev-_b[diffmod]*diffmod-_b[diffage]*diffage-_b[diffage2]*diffage2-_b[diffmarried]*diffmarried-dyreff if e(sample)					

forvalues j=1(1)21	{
	replace yreff_nosel=yreff_nosel+_b[diffyrd`j']*yrd`j'  if e(sample)
					}
	gen reslev_nosel=lw-_b[diffsev]*sev-_b[diffmod]*mod-_b[diffage]*age-_b[diffage2]*age2-_b[diffmarried]*diffmarried-yreff_nosel  if e(sample)


*Overall Wage equation for the household head:
xi: reg difflw diffsev diffmod diffage diffage2   diffmarried diffwhite diffyrd* i.gap|lamda0 i.gap|lag1lamda0 ,vce(cluster newid) nocons
*testparm lamda0-_IgapXlag1_21

replace dyreff=. if !e(sample)
replace yreff=. if !e(sample)

forvalues j=1(1)21	{
	replace dyreff=dyreff+_b[diffyrd`j']*diffyrd`j' if e(sample)
					}
gen g1=difflw-_b[diffsev]*diffsev-_b[diffmod]*diffmod-_b[diffage]*diffage-_b[diffage2]*diffage2-_b[diffmarried]*diffmarried-dyreff 	if e(sample)				

forvalues j=1(1)21	{
	replace yreff=yreff+_b[diffyrd`j']*yrd`j' if e(sample)
					}
	gen reslev=lw-_b[diffsev]*sev-_b[diffmod]*mod-_b[diffage]*age-_b[diffage2]*age2-_b[diffmarried]*diffmarried-yreff_rec  if e(sample)

*Overall wage equation for household head, with 2008+ dummy
xi: reg difflw diffsev diffmod diffage diffage2   diffmarried diffwhite diffyrd* i.gap|lamda0_rec i.gap|lag1lamda0_rec ,vce(cluster newid) nocons
*testparm lamda0-_IgapXlag1_21

replace dyreff_rec=. if !e(sample)
replace yreff_rec=. if !e(sample)

forvalues j=1(1)21	{
	replace dyreff_rec=dyreff_rec+_b[diffyrd`j']*diffyrd`j' if e(sample)
					}
gen g1_rec=difflw-_b[diffsev]*diffsev-_b[diffmod]*diffmod-_b[diffage]*diffage-_b[diffage2]*diffage2-_b[diffmarried]*diffmarried-dyreff_rec 	if e(sample)				

forvalues j=1(1)21	{
	replace yreff_rec=yreff_rec+_b[diffyrd`j']*yrd`j' if e(sample)
					}
	gen reslev_rec=lw-_b[diffsev]*sev-_b[diffmod]*mod-_b[diffage]*age-_b[diffage2]*age2-_b[diffmarried]*diffmarried-yreff_rec  if e(sample)

	
********WAGE EQUATION WITH GROUP FEs IN SELECTION************
*Overall Wage equation for the household head:
*interacting age with severity (vocational grid stuff as a shifter):
g diffsev_age = diffsev*age
g diffmod_age = diffmod*age


g lamda0_fe_age = lamda0_fe*age
g lag1lamda0_fe_age = lag1lamda0_fe*age

g lamda0_fe_dsev = lamda0_fe*diffsev
g lamda0_fe_dmod = lamda0_fe*diffmod
g lamda0_fe_sev = lamda0_fe*sev
g lamda0_fe_mod = lamda0_fe*mod
g lag1lamda0_fe_dsev = lag1lamda0_fe*diffsev
g lag1lamda0_fe_dmod = lag1lamda0_fe*diffmod
g lag1lamda0_fe_sev = lag1lamda0_fe*sev
g lag1lamda0_fe_mod = lag1lamda0_fe*mod

g lamda0_fe_agedsev = lamda0_fe*diffsev*age
g lamda0_fe_agedmod = lamda0_fe*diffmod*age
g lamda0_fe_agesev = lamda0_fe*sev*age
g lamda0_fe_agemod = lamda0_fe*mod*age
g lag1lamda0_fe_agedsev = lag1lamda0_fe*diffsev*age
g lag1lamda0_fe_agedmod = lag1lamda0_fe*diffmod*age
g lag1lamda0_fe_agesev = lag1lamda0_fe*sev*age
g lag1lamda0_fe_agemod = lag1lamda0_fe*mod*age

g sp_lamda0_fe_age = sp_lamda0_fe*sp_age
g lag1sp_lamda0_fe_age = lag1sp_lamda0_fe*sp_age

g sp_lamda0_fe_dmod = sp_lamda0_fe*diffsp_mod
g sp_lamda0_fe_mod = sp_lamda0_fe*sp_mod
g lag1sp_lamda0_fe_dmod = lag1sp_lamda0_fe*diffsp_mod
g lag1sp_lamda0_fe_mod = lag1sp_lamda0_fe*sp_mod

g sp_lamda0_fe_agedmod = sp_lamda0_fe*diffsp_mod*sp_age
g sp_lamda0_fe_agemod = sp_lamda0_fe*sp_mod*sp_age
g lag1sp_lamda0_fe_agedmod = lag1sp_lamda0_fe*diffsp_mod*sp_age
g lag1sp_lamda0_fe_agemod = lag1sp_lamda0_fe*sp_mod*sp_age


*Simpler specification:
xi: reg difflw diffsev diffmod diffage diffage2   diffmarried diffwhite diffyrd* i.gap|lamda0_fe i.gap|lag1lamda0_fe ,vce(cluster newid) nocons

reg lw prodtype_* sev mod age age2 old sev_old mod_old married yrd* , vce(cluster newid) nocons

predict lw_res, residuals

eststo headresults_model

xi: reg difflw diffsev diffmod diffage diffage2   diffmarried diffwhite diffyrd* i.gap|lamda0_fe i.gap|lag1lamda0_fe i.gap|lamda0_fe_age i.gap|lag1lamda0_fe_age i.gap|lamda0_fe_sev i.gap|lag1lamda0_fe_sev i.gap|lamda0_fe_mod i.gap|lag1lamda0_fe_mod i.gap|lamda0_fe_dsev i.gap|lag1lamda0_fe_dsev i.gap|lamda0_fe_dmod i.gap|lag1lamda0_fe_dmod i.gap|lamda0_fe_agesev i.gap|lag1lamda0_fe_agesev i.gap|lamda0_fe_agemod i.gap|lag1lamda0_fe_agemod i.gap|lamda0_fe_agedsev i.gap|lag1lamda0_fe_agedsev i.gap|lamda0_fe_agedmod i.gap|lag1lamda0_fe_agedmod,vce(cluster newid) nocons

eststo headresults


testparm lamda0_fe-_IgapXlag1_2
*test _b[lamda0_fe] = _b[_IgapXlamda_2]
replace dyreff_fe=. if !e(sample)

forvalues j=1(1)21	{
	replace dyreff_fe=dyreff_fe+_b[diffyrd`j']*diffyrd`j' if e(sample)
					}
gen g1_fe=difflw-_b[diffsev]*diffsev-_b[diffmod]*diffmod-_b[diffage]*diffage-_b[diffage2]*diffage2-_b[diffmarried]*diffmarried-dyreff_fe 	if e(sample)				
forvalues j=1(1)21	{
	replace yreff_fe=yreff_fe+_b[diffyrd`j']*yrd`j' if e(sample)
					}
	gen reslev_fe=lw-_b[diffsev]*sev-_b[diffmod]*mod-_b[diffage]*age-_b[diffage2]*age2-_b[diffmarried]*diffmarried-yreff_fe if e(sample)
*removing selection term from the individual residuals:
*For correcting the FE, I use non-lagged lamda0_fe b/c that represents selection in today's wages from choosing to work today.
*Only correcting for selection on current wage shocks, not past ones.
gen correct_fe = - _b[lamda0_fe]*lamda0_fe

*gen correct_fe = - ((age-22)^0.5)*_b[lamda0_fe]*lamda0_fe
*Why taking square root of age?
*gen correct_fe = - ((age-22))*_b[lamda0_fe]*lamda0_fe
*ALTERNATIVE: use parameters estimated below
*gen correct_fe = - ((age-22))*lamda0_fe*(0.0254078^0.5)*(0.7455936)
*gen correct_fe = - ((age-22-1))*lamda0_fe*(0.0254078^0.5)*(-0.8849552) - lamda0_fe*(0.0254078^0.5)*(0.7455936)

replace reslev_fe = reslev_fe + correct_fe

*Separately by marital status
forvalues m = 0(1)1{
	xi: reg difflw diffsev diffmod diffage diffage2   diffmarried diffwhite diffyrd* i.gap|lamda0_fem`m' i.gap|lag1lamda0_fem`m' ,vce(cluster newid) nocons

}

xi: reg lw prodtype_*, vce(cluster newid) nocons
eststo headmeans

*Wage equation for spouses:
label var diffsp_sev "Health - Severe"
label var diffsp_mod "Health - Moderate"
label var diffsp_age "Age"
label var diffsp_age2 "Age Squared"

set scheme s1mono

*without selection correction:
xi: reg diffsp_lw diffsp_dis diffsp_age diffsp_age2    diffwhite  diffyrd* if married==1 ,vce(cluster newid) nocons
eststo spouseresults_noselection


replace sp_dyreff_nosel=. if !e(sample)
replace sp_g1_nosel = . if !e(sample)
replace sp_yreff_nosel=. if !e(sample)
replace sp_reslev_nosel=. if !e(sample)

forvalues j=1(1)21	{
	replace sp_dyreff_nosel=sp_dyreff_nosel+_b[diffyrd`j']*diffyrd`j' if e(sample)
					}
replace sp_g1_nosel=diffsp_lw-_b[diffsp_dis]*diffsp_dis-_b[diffsp_age]*diffsp_age-_b[diffsp_age2]*diffsp_age2-sp_dyreff if e(sample)						

forvalues j=1(1)21	{
	replace sp_yreff_nosel=sp_yreff_nosel+_b[diffyrd`j']*yrd`j' if e(sample)
					}
	replace sp_reslev_nosel=lw-_b[diffsp_dis]*sp_dis-_b[diffsp_age]*sp_age-_b[diffsp_age2]*sp_age2-sp_yreff_nosel if e(sample)


*with selection correction:
xi: reg diffsp_lw diffsp_dis diffsp_age diffsp_age2    diffwhite  diffyrd* i.gap|sp_lamda0 i.gap|lag1sp_lamda0 if married==1 ,vce(cluster newid) nocons
coefplot , keep(diffsp_dis diffsp_age diffsp_age2) xline(0, lcolor(red) lpattern(dash))
*graph export "wagecoefs_spouse.pdf", replace

gen sp_emp = sp_lw != .
gen emp = lw != .

gen sp_diffemp = diffsp_lw != .
gen diffemp = difflw != .

replace sp_dyreff=. if !e(sample)
replace sp_g1 = . if !e(sample)
replace sp_yreff=. if !e(sample)
replace sp_reslev=. if !e(sample)

replace sp_dyreff = sp_dyreff 
forvalues j=1(1)21	{
	replace sp_dyreff=sp_dyreff+_b[diffyrd`j']*diffyrd`j' if e(sample)
					}
replace sp_g1=diffsp_lw-_b[diffsp_dis]*diffsp_dis-_b[diffsp_age]*diffsp_age-_b[diffsp_age2]*diffsp_age2-sp_dyreff if e(sample)

forvalues j=1(1)21	{
	replace sp_yreff=sp_yreff+_b[diffyrd`j']*yrd`j' if e(sample)
					}
	replace sp_reslev=sp_lw-_b[diffsp_dis]*sp_dis-_b[diffsp_age]*sp_age-_b[diffsp_age2]*sp_age2-sp_yreff if e(sample)
	
	*************SELECTION CORRECTION AND GROUP FE********

	
reg sp_lw prodtype_*  sp_mod sp_age sp_age2 sp_old sp_mod_old yrd1-yrd21, vce(cluster newid) nocons
eststo spouseresults_model

predict sp_lw_res, residuals

sort newid year 
qby newid: g laglw_res=lw_res[_n-1]
qby newid: g lagsp_lw_res=lw_res[_n-1]

sum lw_res
estadd scalar wagevar = `r(sd)'^2

sum sp_lw_res
estadd scalar spwagevar = `r(sd)'^2

corr lw_res laglw_res, covariance
estadd scalar wagecov = `r(cov_12)' 

corr sp_lw_res lagsp_lw_res, covariance
estadd scalar spwagecov = `r(cov_12)'

corr lw_res sp_lw_res
estadd scalar wagecorr = `r(rho)'


gen lw_res2 = lw_res*lw_res
gen sp_lw_res2 = sp_lw_res*sp_lw_res
gen lwcov = lw_res*sp_lw_res

foreach x in "" "sp_" {
gen `x'lwlaglw_res = `x'lw_res*lag`x'lw_res
sum `x'lw_res if lag`x'lw_res != .
local mean1 = `r(mean)'
sum lag`x'lw_res if `x'lw_res != .
local mean2 = `r(mean)'
replace `x'lwlaglw_res = `x'lwlaglw_res - `mean1'*`mean2'
}

*Getting wage variances and autocovariances
gmm (eq1: {wagevar=1.0} - lw_res2) ///
	(eq2: {wagecov=1.0} - lwlaglw_res) ///
	(eq3: {spwagevar=1.0} - sp_lw_res2) ///
	(eq4: {spwagecov=1.0} - sp_lwlaglw_res), nocommonesample winitial(identity) vce(cluster newid) wmat(cluster newid)

	eststo wagevars

			esttab wagevars using "wagevars`1'_jump.tex", keep(wagevar:_cons wagecov:_cons spwagevar:_cons spwagecov:_cons) coeflabels() noeqlines eqlabels("Variance - Head" "Autocovariance - Head" "Variance - Spouse" "Autocovariance - Spouse", none) se noobs substitute(\_ _) replace nonum nolines nonotes fragment nomtitles
		estout wagevars using "wagevars`1'_jump.txt", keep(wagevar:_cons wagecov:_cons spwagevar:_cons spwagecov:_cons)  replace nonum eqlabels(, none)  cells("b se p")  mlabels(,none)  

*Getting wagecorr (note common sample):
gmm (eq1: {wagevar=1.0} - lw_res2) ///
	(eq2: {spwagevar=1.0} - sp_lw_res2) ///
	(eq3: {wagecorr=0.0} - lwcov), nocommonesample winitial(identity) vce(cluster newid) wmat(cluster newid)

		eststo wagecorrs

				esttab wagecorrs using "wagevars`1'_jump.tex", keep(wagecorr:_cons) coeflabels() noeqlines eqlabels("intracorr", none) varlabels(_cons wagecorr) se noobs substitute(\_ _) append nonum nolines nonotes fragment nomtitles
		estout wagecorrs using "wagevars`1'_jump.txt", keep(wagecorr:_cons)  append nonum eqlabels("intracorr",none) varlabels(_cons intracorr) collabels(, none) cells("b se p")  mlabels(,none)  

*with selection correction:
preserve 
drop diffmod diffage diffage2 mod age age2 old mod_old
ren diffsp_dis diffmod
ren diffsp_age diffage
ren diffsp_age2 diffage2
ren sp_age age
ren sp_age2 age2
ren sp_mod mod
ren sp_old old
ren sp_mod_old mod_old

xi: reg diffsp_lw diffmod diffage diffage2    diffwhite  diffyrd* i.gap|sp_lamda0_fe i.gap|lag1sp_lamda0_fe if married==1,vce(cluster newid) nocons

xi: reg diffsp_lw diffmod diffage diffage2   diffwhite diffyrd* i.gap|sp_lamda0_fe i.gap|lag1sp_lamda0_fe i.gap|sp_lamda0_fe_age i.gap|lag1sp_lamda0_fe_age i.gap|sp_lamda0_fe_mod i.gap|lag1sp_lamda0_fe_mod i.gap|sp_lamda0_fe_dmod i.gap|lag1sp_lamda0_fe_dmod i.gap|sp_lamda0_fe_agemod i.gap|lag1sp_lamda0_fe_agemod i.gap|sp_lamda0_fe_agedmod i.gap|lag1sp_lamda0_fe_agedmod,vce(cluster newid) nocons

testparm sp_lamda0_fe-_IgapXlag1_2

eststo spouseresults

xi: reg diffsp_lw diffmod diffage diffage2   diffwhite diffyrd* i.gap|sp_lamda0_fe i.gap|lag1sp_lamda0_fe i.gap|sp_lamda0_fe_age i.gap|lag1sp_lamda0_fe_age i.gap|sp_lamda0_fe_mod i.gap|lag1sp_lamda0_fe_mod i.gap|sp_lamda0_fe_dmod i.gap|lag1sp_lamda0_fe_dmod i.gap|sp_lamda0_fe_agemod i.gap|lag1sp_lamda0_fe_agemod i.gap|sp_lamda0_fe_agedmod i.gap|lag1sp_lamda0_fe_agedmod,vce(cluster newid) nocons



xi: reg sp_lw prodtype_*, vce(cluster newid) nocons
eststo spousemeans

restore
xi: reg diffsp_lw diffsp_dis diffsp_age diffsp_age2    diffwhite  diffyrd* i.gap|sp_lamda0_fe i.gap|lag1sp_lamda0_fe if married==1,vce(cluster newid) nocons

xi: reg diffsp_lw diffsp_dis diffsp_age diffsp_age2   diffwhite diffyrd* i.gap|sp_lamda0_fe i.gap|lag1sp_lamda0_fe i.gap|sp_lamda0_fe_age i.gap|lag1sp_lamda0_fe_age i.gap|sp_lamda0_fe_mod i.gap|lag1sp_lamda0_fe_mod i.gap|sp_lamda0_fe_dmod i.gap|lag1sp_lamda0_fe_dmod i.gap|sp_lamda0_fe_agemod i.gap|lag1sp_lamda0_fe_agemod i.gap|sp_lamda0_fe_agedmod i.gap|lag1sp_lamda0_fe_agedmod,vce(cluster newid) nocons

coefplot , keep(diffsp_dis diffsp_age diffsp_age2) xline(0, lcolor(red) lpattern(dash))
graph export "wagecoefs_spouse`1'_jump.pdf", replace




replace sp_dyreff_fe=. if !e(sample)
replace sp_g1_fe = . if !e(sample)
replace sp_yreff_fe = . if !e(sample)
replace sp_reslev_fe = . if !e(sample)

replace sp_dyreff_fe = sp_dyreff_fe 
forvalues j=1(1)21	{
	replace sp_dyreff_fe=sp_dyreff_fe+_b[diffyrd`j']*diffyrd`j' if e(sample)
					}
replace sp_g1_fe=diffsp_lw-_b[diffsp_dis]*diffsp_dis-_b[diffsp_age]*diffsp_age-_b[diffsp_age2]*diffsp_age2-sp_dyreff_fe if e(sample)
	
forvalues j=1(1)21	{
	replace sp_yreff_fe=sp_yreff_fe+_b[diffyrd`j']*yrd`j' if e(sample)
					}
replace sp_reslev_fe=sp_lw-_b[diffsp_dis]*sp_dis-_b[diffsp_age]*sp_age-_b[diffsp_age2]*sp_age2-sp_yreff_fe if e(sample)
*removing selection term from the individual residuals:
*For correcting the FE, I use non-lagged lamda0_fe b/c that represents selection in today's wages from choosing to work today.
*Only correcting for selection on current wage shocks, not past ones.
gen sp_correct_fe = - _b[sp_lamda0_fe]*sp_lamda0_fe

*gen sp_correct_fe = - ((age-22)^0.5)*_b[sp_lamda0_fe]*sp_lamda0_fe
*Why taking square root of age?
*gen sp_correct_fe = - ((age-22))*_b[sp_lamda0_fe]*sp_lamda0_fe
*ALTERNATIVE: use parameters estimated below
*gen sp_correct_fe = - ((age-22))*sp_lamda0_fe*(sp_sigmaz^0.5)*(rho_00)
*gen sp_correct_fe = - ((age-22-1))*sp_lamda0_fe*(sp_sigmaz^0.5)*(-rho_01) - sp_lamda0_fe*(sp_sigmaz^0.5)*(rho_00)

replace sp_reslev_fe = sp_reslev_fe + sp_correct_fe

****************************************************************************************

egen fix=mean(reslev),by(newid)
egen fix_fe = mean(reslev_fe), by(type_prodwage)
egen fix_nosel=mean(reslev_nosel),by(newid)
egen sp_fix=mean(sp_reslev),by(sp_newid)
egen sp_fix_nosel=mean(sp_reslev_nosel),by(sp_newid)
egen sp_fix_fe = mean(sp_reslev_fe), by(type_prodwage)

egen sp_intercept = mean(sp_yreff)
egen intercept = mean(yreff)

egen sp_intercept_fe = mean(sp_yreff_fe)
egen intercept_fe = mean(yreff_fe)

g g2=g1 if gap==2
g g2_fe=g1_fe if gap==2
replace g1=. if gap!=1
replace g1_fe=. if gap!=1

g g2_nosel=g1_nosel if gap==2
replace g1_nosel=. if gap!=1

g sp_g2=sp_g1 if gap==2
replace sp_g1=. if gap!=1

g sp_g2_fe=sp_g1_fe if gap==2
replace sp_g1_fe=. if gap!=1

g sp_g2_nosel=sp_g1_nosel if gap==2
replace sp_g1_nosel=. if gap!=1

bysort type_prodwage: sum fix_fe sp_fix_fe
bysort type_prodwage: sum sp_reslev_fe reslev_fe intercept_fe sp_intercept_fe  correct_fe sp_correct_fe
preserve
duplicates drop newid, force
bysort type_prodwage: sum sp_reslev_fe reslev_fe intercept_fe sp_intercept_fe  correct_fe sp_correct_fe
restore


sum intercept_fe
local constant = `r(mean)'
sum sp_intercept_fe
local sp_constant = `r(mean)'
levelsof type_prodwage, local(types)
*printing regression results:
local stat_types
foreach i of local types{
estimates restore headresults
sum reslev_fe if type_prodwage == `i'
estadd scalar Type`i' = `r(mean)' + `constant'

estimates restore spouseresults
sum sp_reslev_fe if type_prodwage == `i'
estadd scalar Type`i' = `r(mean)' + `sp_constant'

local stat_types `stat_types' Type`i'
}
esttab headresults spouseresults using "wageparams`1'_jump.tex", tex fragment keep(diffage diffage2 diffmod diffsev diffmarried) scalars(Type1 Type2 Type3) coeflabels(diffsev "Sev. Dis." diffmod "Mod. Dis." diffage "Age " diffage2 " $ Age^2/100$" diffmarried "Married") replace mtitles("Head" "Partner") nonum se stats(Type1 Type2 Type3, labels("$ f_1$" "$ f_2$" "$ f_3$")) substitute(\_ _)
estout headresults using "wageparams`1'_jump.txt", keep(diffage diffage2 diffmod diffsev diffmarried)  replace nonum eqlabels(, none)  cells("b se p")  mlabels(,none)  stats(`stat_types')
estout headresults_model using "wagemoments`1'_jump.txt", replace nonum eqlabels(, none)  cells("b se p")  mlabels(,none)  
*estout headmeans using "wagemeans`1'_jump.txt", replace nonum eqlabels(, none)  cells("b se p")  mlabels(,none)  
estout headresults_noselection using "wageparams_noselection`1'.txt", keep(diffage diffage2 diffmod diffsev diffmarried)  replace nonum eqlabels(, none)  cells("b se p")  mlabels(,none)  stats(`stat_types')
estout spouseresults using "spwageparams`1'_jump.txt", keep(diffage diffage2 diffmod)  replace nonum eqlabels(, none)  cells("b se p")  mlabels(,none)  stats(`stat_types')
estout spouseresults_noselection using "spwageparams_noselection`1'.txt", keep(diffsp_age diffsp_age2 diffsp_dis)  replace nonum eqlabels(, none)  cells("b se p")  mlabels(,none)
estout spouseresults_model using "spwagemoments`1'_jump.txt", replace nonum eqlabels(, none)  cells("b se p")  mlabels(,none)  
*estout spouseresults_model using "wagevars`1'_jump.txt", drop(*) replace nonum eqlabels(, none)  cells("b se p")  mlabels(,none)  stats(wagevar wagecov spwagevar spwagecov wagecorr)
*estout spousemeans using "spwagemeans`1'_jump.txt", replace nonum eqlabels(, none)  cells("b se p")  mlabels(,none)  


tempfile sp_temp_fix
save `sp_temp_fix',replace


gen g1sq=g1^2
gen g2sq=g2^2

gen g1sq_fe=g1_fe^2
gen g2sq_fe=g2_fe^2

gen g1sq_nosel=g1_nosel^2
gen g2sq_nosel=g2_nosel^2

gen sp_g1sq=sp_g1^2
gen sp_g2sq=sp_g2^2

gen sp_g1sq_fe=sp_g1_fe^2
gen sp_g2sq_fe=sp_g2_fe^2

gen sp_g1sq_nosel=sp_g1_nosel^2
gen sp_g2sq_nosel=sp_g2_nosel^2

ren lamda0 l_p
ren lamda0_fe l_p_fe
ren both_lamda0_fe both_l_p_fe
ren lag1lamda0 lagl_p
ren lag1lamda0_fe lagl_p_fe
ren both_lag1lamda0_fe both_lagl_p_fe
ren xb0 z_g
ren xb0_fe z_g_fe
ren both_xb0_fe both_z_g_fe
ren lag1xb0 lagz_g
ren lag1xb0_fe lagz_g_fe
ren both_lag1xb0_fe both_lagz_g_fe

ren sp_lamda0 sp_l_p
ren sp_lamda0_fe sp_l_p_fe
ren lag1sp_lamda0 laglsp_p
ren lag1sp_lamda0_fe laglsp_p_fe
ren sp_xb0 sp_z_g
ren sp_xb0_fe sp_z_g_fe
ren lag1sp_xb0 lagsp_z_g
ren lag1sp_xb0_fe lagsp_z_g_fe


keep newid married year *g1* *g2* *l_p* lagl*_p* *z_g* lag*z_g* age sev mod gap sp_sev sp_mod sp_age birthy sp_birthy
preserve
keep if gap==1|gap==2
gen g=g1 if gap==1
replace g=g2 if gap==2

gen g_fe = g1_fe if gap==1
replace g_fe = g2_fe if gap==2

gen g1_both_fe = g1_fe if gap == 1 & sp_g1_fe != .
gen sp_g1_both_fe = sp_g1_fe if gap == 1 & g1_fe != .

gen gsp_g1_fe = g1_fe*sp_g1_fe 

drop g1 g2 
drop g1_fe g2_fe

gen g2=g1sq if gap==1
replace g2=g2sq if gap==2

gen g2_fe = g1sq_fe if gap==1
replace g2_fe = g2sq_fe if gap==2

gen sp_g=sp_g1 if gap==1
replace sp_g=sp_g2 if gap==2

gen sp_g_fe=sp_g1_fe if gap==1
replace sp_g_fe=sp_g2_fe if gap==2

gen gsp_g_fe = g_fe*sp_g_fe 

drop sp_g1 sp_g2
drop sp_g1_fe sp_g2_fe

gen sp_g2=sp_g1sq if gap==1
replace sp_g2=sp_g2sq if gap==2

gen sp_g2_fe=sp_g1sq_fe if gap==1
replace sp_g2_fe=sp_g2sq_fe if gap==2

g gap1=gap==1
g gap2=gap==2
keep newid married year both* l_p* lagl_p* z_g* lagz_g* age sev mod gap1 gap2 g* g2* birthy sp_birthy sp_l_p* laglsp_p* sp_z_g* lagsp_z_g* sp_age sp_sev sp_mod sp_g* sp_g2* 
sort newid year 
qui by newid:gen gg1=(g)*(g[_n-1])
sort newid year 
qui by newid:gen gg2=(g)*(g[_n-2])

qui by newid:gen gg1_fe=(g_fe)*(g_fe[_n-1])

g lpg1=l_p*gap1
g lpg1_fe=l_p_fe*gap1
g lpg1both_fe = both_l_p_fe*gap1
g lpg2both_fe = both_l_p_fe*gap2
g laglpg1=lagl_p*gap1
g laglpg1_fe=lagl_p_fe*gap1
g laglpg1both_fe=both_lagl_p_fe*gap1
g laglpg2both_fe=both_lagl_p_fe*gap2
g lpg2=l_p*gap2
g lpg2_fe=l_p_fe*gap2
g laglpg2=lagl_p*gap2
g laglpg2_fe=lagl_p_fe*gap2

g zglpg1=z_g*l_p*gap1
g zglpg1_fe=z_g_fe*l_p_fe*gap1
g zglpg1both_fe=both_z_g_fe*both_l_p_fe*gap1
g lagzglaglpg1=lagz_g*lagl_p*gap1
g lagzglaglpg1_fe=lagz_g_fe*lagl_p_fe*gap1
g lagzglaglpg1both_fe=both_lagz_g_fe*both_lagl_p_fe*gap1
g zglpg2=z_g*l_p*gap2
g zglpg2_fe=z_g_fe*l_p_fe*gap2
g lagzglaglpg2=lagz_g*lagl_p*gap2
g lagzglaglpg2_fe=lagz_g_fe*lagl_p_fe*gap2

sort newid year 
qui by newid:gen sp_gg1=(sp_g)*(sp_g[_n-1])
sort newid year 
qui by newid:gen sp_gg2=(sp_g)*(sp_g[_n-2])

qui by newid: gen sp_gg1_fe=(sp_g_fe)*(sp_g_fe[_n-1])

g lsp_pg1=sp_l_p*gap1
g lsp_pg1_fe=sp_l_p_fe*gap1
g laglsp_pg1=laglsp_p*gap1
g laglsp_pg1_fe=laglsp_p_fe*gap1
g lsp_pg2=sp_l_p*gap2
g lsp_pg2_fe=sp_l_p_fe*gap2
g laglsp_pg2=laglsp_p*gap2
g laglsp_pg2_fe=laglsp_p_fe*gap2

g sp_zglpg1=sp_z_g*sp_l_p*gap1
g sp_zglpg1_fe=sp_z_g_fe*sp_l_p_fe*gap1
g lagsp_zglaglsp_pg1=lagsp_z_g*laglsp_p*gap1
g lagsp_zglaglsp_pg1_fe=lagsp_z_g_fe*laglsp_p_fe*gap1
g sp_zglpg2=sp_z_g*sp_l_p*gap2
g sp_zglpg2_fe=sp_z_g_fe*sp_l_p_fe*gap2
g lagsp_zglaglsp_pg2=lagsp_z_g*laglsp_p*gap2
g lagsp_zglaglsp_pg2_fe=lagsp_z_g_fe*laglsp_p_fe*gap2

****ADDRESS VARIANCE HERE*****

		//MY VERSION OF ESTIMATED VARIANCE, MY WAGE EQUATION:
		///*AND FIXING THE SELECTION CORRECTION IN EQ2:
		*rho_ts_ij is corr b/w xi of partner t and nu of partner s with j lags on nu or i lags on xi
		//mu1 and mu2 are nuisance parameters, non-zero mean of g_fe and sp_g_fe
		//theta is correlation term

gmm  ///
		(eq1: g_fe - {mu1} - sqrt({sigma_z1=0.2})*{rho_11_00=0.0}*lpg1_fe - sqrt({sigma_z1})*{rho_11_01=0.0}*laglpg1_fe	///
				- sqrt({sigma_z1})*({rho_11_00}+{rho_11_10=0.0})*lpg2_fe - sqrt({sigma_z1})*({rho_11_01}+{rho_11_02=0.0})*laglpg2_fe) ///
		(eq2: g2_fe- ({sigma_z1}+2*{sigma_e1=0.1})*gap1 ///
		           - (2*{sigma_z1}+2*{sigma_e1})*gap2 ///
				+ {sigma_z1}*{rho_11_00}^2*(zglpg1_fe)+ {sigma_z1}*{rho_11_01}^2*(lagzglaglpg1_fe)	///
				+ {sigma_z1}*({rho_11_00}^2+{rho_11_10}^2)*(zglpg2_fe) + {sigma_z1}*({rho_11_01}^2+{rho_11_02}^2)*(lagzglaglpg2_fe) ///
				- 2*{sigma_z1}*lpg1_fe*laglpg1_fe*{rho_11_00}*{rho_11_01} ///
				- 2*{sigma_z1}*lpg2_fe*laglpg2_fe*({rho_11_00}*{rho_11_02}+{rho_11_10}*{rho_11_01}) ///
				- 2*{sigma_z1}*((lpg2_fe*{rho_11_00}+laglpg2_fe*{rho_11_02})*(lpg2_fe*{rho_11_10}+laglpg2_fe*{rho_11_01})) 	///
			) ///
		(eq3: gg1_fe + {sigma_e1})	///
		(eq4: sp_g_fe - {mu2} - sqrt({sigma_z2=0.2})*{rho_22_00=0.0}*lsp_pg1_fe - sqrt({sigma_z2})*{rho_22_01=0.0}*laglsp_pg1_fe	///
				- sqrt({sigma_z2})*({rho_22_00}+{rho_22_10=0.0})*lsp_pg2_fe - sqrt({sigma_z2})*({rho_22_01}+{rho_22_02=0.0})*laglsp_pg2_fe) ///
		(eq5: sp_g2_fe- ({sigma_z2}+2*{sigma_e2=0.1})*gap1 -(2*{sigma_z2}+2*{sigma_e2})*gap2 ///
				+ {sigma_z2}*{rho_22_00}^2*(sp_zglpg1_fe)+ {sigma_z2}*{rho_22_01}^2*(lagsp_zglaglsp_pg1_fe)	///
				+ {sigma_z2}*({rho_22_00}^2+{rho_22_10}^2)*(sp_zglpg2_fe) + {sigma_z2}*({rho_22_01}^2+{rho_22_02}^2)*(lagsp_zglaglsp_pg2_fe) ///
				- 2*{sigma_z2}*lsp_pg1_fe*laglsp_pg1_fe*{rho_22_00}*{rho_22_01} ///
				- 2*{sigma_z2}*lsp_pg2_fe*laglsp_pg2_fe*({rho_22_00}*{rho_22_02}+{rho_22_10}*{rho_22_01}) ///
				- 2*{sigma_z2}*((lsp_pg2_fe*{rho_22_00}+laglsp_pg2_fe*{rho_22_02})*(lsp_pg2_fe*{rho_22_10}+laglsp_pg2_fe*{rho_22_01})) 	///				
				) ///
		(eq6: sp_gg1_fe + {sigma_e2})	///
		(eq7: gsp_g_fe - sqrt({sigma_z1})*sqrt({sigma_z2})*( /// INTRA-HOUSEHOLD CORRELATION IN SHOCKS, WITH SELECTION
					 ({rho_11_00}*lpg1_fe+{rho_11_01}*laglpg1_fe)*({rho_22_00}*lsp_pg1_fe+{rho_22_01}*laglsp_pg1_fe) /// 
 				   + {theta_z}*gap1) /// 
				   - sqrt({sigma_z1})*sqrt({sigma_z2})*( /// TWO YEAR DIFF:
				   (({rho_11_00}+{rho_11_10})*lpg2_fe + ({rho_11_02}+{rho_11_01})*laglpg2_fe)*  ///
				   (({rho_22_00}+{rho_22_10})*lsp_pg2_fe +({rho_22_02}+{rho_22_01})*laglsp_pg2_fe) ///
				   + 2*{theta_z}*gap2) ///
						- 2*sqrt({sigma_e1})*sqrt({sigma_e2})*{theta_e} ///
		) ///
		, instruments (eq1: lpg1_fe laglpg1_fe lpg2_fe laglpg2_fe) ///
		instruments (eq2: gap1 gap2 zglpg1_fe lagzglaglpg1_fe zglpg2_fe lagzglaglpg2_fe) ///
		instruments (eq4: lsp_pg1_fe laglsp_pg1_fe lsp_pg2_fe laglsp_pg2_fe) ///
		instruments (eq5: gap1 gap2 sp_zglpg1_fe lagsp_zglaglsp_pg1_fe sp_zglpg2_fe lagsp_zglaglsp_pg2_fe) ///
		instruments (eq7: lpg1both_fe laglpg1both_fe lpg2both_fe laglpg2both_fe) ///
		nocommonesample winitial(identity) vce(cluster newid) wmat(cluster newid)
		eststo varcov

		esttab varcov using "wagevarcov`1'_jump.tex", keep(sigma_z1:_cons sigma_z2:_cons theta_z:_cons sigma_e1:_cons sigma_e2:_cons) coeflabels() noeqlines eqlabels(\$\sigma_z^1\$ \$\sigma_z^2\$ \$\rho_z\$ \$\sigma_e^1\$ \$\sigma_e^2\$, none) se noobs substitute(\_ _) replace nonum nolines nonotes fragment nomtitles
		estout varcov using "wagevarcov`1'_jump.txt", keep(sigma_z1:_cons sigma_z2:_cons theta_z:_cons)  replace nonum eqlabels(, none)  cells("b se p")  mlabels(,none)  

restore 

preserve


*head variance, no selection:
drop g1 
drop g2
rename g1_nosel g1
rename g2_nosel g2
keep if gap==1|gap==2
gen g=g1 if gap==1
replace g=g2 if gap==2
drop g1 g2
gen g2=g1sq_nosel if gap==1
replace g2=g2sq_nosel if gap==2
g gap1=gap==1
g gap2=gap==2
keep newid year age sp_sev mod gap1 gap2 g g2 birthy sp_birthy
sort newid year 
qui by newid:gen gg1=(g)*(g[_n-1])
sort newid year 
qui by newid:gen gg2=(g)*(g[_n-2])

gmm 	(eq2: g2- ({sigma_z}+2*{sigma_e=0.25})*gap1 -(2*{sigma_z}+2*{sigma_e})*gap2)	///
		(eq3: gg1 + {sigma_e})	, ///
			 nocommonesample winitial(identity) vce(cluster newid) wmat(cluster newid)





restore
preserve
*spousal wage variance, no selection:
drop sp_g1 
drop sp_g2
rename sp_g1_nosel sp_g1
rename sp_g2_nosel sp_g2
keep if gap==1|gap==2
gen sp_g=sp_g1 if gap==1
replace sp_g=sp_g2 if gap==2
drop sp_g1 sp_g2
gen sp_g2=sp_g1sq_nosel if gap==1
replace sp_g2=sp_g2sq_nosel if gap==2
g gap1=gap==1
g gap2=gap==2
keep newid year sp_l_p laglsp_p sp_z_g lagsp_z_g sp_age sp_sev mod gap1 gap2 sp_g sp_g2 birthy sp_birthy
sort newid year 
qui by newid:gen sp_gg1=(sp_g)*(sp_g[_n-1])
sort newid year 
qui by newid:gen sp_gg2=(sp_g)*(sp_g[_n-2])

gmm (eq2: sp_g2- ({sigma_z}+2*{sigma_e=0.25})*gap1 -(2*{sigma_z}+2*{sigma_e})*gap2)	///
		(eq3: sp_gg1 + {sigma_e})	, ///
			 nocommonesample winitial(identity) vce(cluster newid) wmat(cluster newid)


restore
clear




****Unobserved productivity Type *****************

**************************************************

u `sp_temp_fix',clear
drop married
bysort newid: egen temp = max(marit==1) if age >= 30
 bysort newid: egen married=max(temp)
 drop temp
 replace married=0 if married==.
 preserve
collapse fix fix_nosel fix_fe intercept intercept_fe ,by(married newid educ birthy  sp_todrop)
su fix  
xtile type_fix=fix ,nq(4)
replace type_fix=2 if type_fix >= 2 & type_fix <=3
replace type_fix=3 if type_fix ==4
replace type_fix=1 if fix==.						/*This assumes that people who are never observed at work comne from the bottom of the f distribution*/
sum fix if type_fix==1
replace fix = r(mean) if fix==. 

bysort type_fix: sum fix intercept 
bysort type_fix: sum fix intercept if married==1 


tempfile headfe
save `headfe', replace

keep newid type_fix married
sort newid

tempfile sp_temp_type
save `sp_temp_type',replace
restore
collapse sp_fix sp_fix_nosel sp_intercept, by(newid sp_newid married sp_educ sp_birthy sp_todrop)
merge m:1 newid using `headfe'

twoway scatter fix sp_fix  || function y = x, range(0 4)

su sp_fix if married==1 
*replace sp_fix=r(min) if sp_fix==. & married==1					
xtile sp_type_fix=sp_fix ,nq(4)
replace sp_type_fix=2 if (sp_type_fix>=2 & sp_type_fix<=3)
replace sp_type_fix=3 if (sp_type_fix==4)

/*
levelsof type_fix, local(types)
foreach x of local types {
sum sp_fix if type_fix==`x' & married==1
replace sp_fix = r(mean) if type_fix==`x' & sp_fix==. & married==1
}
*/
keep newid sp_newid sp_type_fix* type_fix* fix sp_fix *intercept married sp_birthy sp_todrop
sort newid

gen sp_fix_interp = sp_fix==.
sum sp_fix, det
replace sp_fix = r(p10) if sp_fix == . & married==1 & sp_todrop ==0
preserve
keep newid sp_newid sp_type_fix married
tempfile sp_temp_type_spouse
save `sp_temp_type_spouse',replace
*type-specific fixed effects:
restore
bysort type_fix: sum sp_fix sp_intercept if married==1 

*40% of married heads never have a spouse with an identified FE.
**************************************************
replace sp_fix = -2 if married==0
quietly sum fix if married==1, det
twoway (scatter  sp_fix fix if married==1, mcolor(dknavy) msize(0.7))  ( scatter sp_fix fix if married==0, mcolor(maroon) jitter(5) msize(0.5)), xline(`r(p50)', lcolor(red) lpattern(dash)) legend(label(1 "Married") label(2 "Unmarried")) xtitle("Head Fixed Effect") ytitle("Spouse Fixed Effect")
graph export "type_corrs_jump.pdf", replace
clear


local type type_model
cd "$revdta"
use  "Full_Panel_2b.dta", clear
merge m:1 newid using "./intermediary/sp_temp_types"

keep if sample_analytic==1


sort newid year
qby newid:gen grly=(exp(sp_lw)/exp(sp_lw[_n-1]))-1
replace grly=. if grly==-1
gen suspect=0
qui su grly,d

replace suspect=grly<-lowlim|(grly>highlim & grly!=.)	
egen ss=sum(suspect),by(newid)				/*# of outlier wage growth records*/
egen sw=sum(exp(sp_lw)<0.5*minwage),by(newid)	/*#of wages below state-level minimum wage*/
egen numy=sum(newid!=.),by(newid)			/*# of obs*/

gen sp_todrop=sw>0|numy<3|ss>0				/*Drop those with outlier wage growth, below 1/2 of mw, or with less than 3 obs*/

cd "$out"

	gen fullearn = exp(lw)*2000 
	gen sp_fullearn = exp(sp_lw)*2000 if sp_todrop==0
	collapse (median) fullearn sp_fullearn, by(age `type')
twoway (line fullearn age if `type'==1, lcolor(black)) (line fullearn age if `type'==2, lcolor(black)) (line fullearn age if `type'==3, lcolor(black)) ///
(line sp_fullearn age if `type'==1, lcolor(red)) (line sp_fullearn age if `type'==2, lcolor(red)) (line sp_fullearn age if `type'==3, lcolor(red)), legend(off) xtitle(Age, size(large)) xlabel(,labsize(large)) ytitle("Earnings (full-time work)", size(large)) ylabel(0(20000)60000,labsize(large)) yscale(range(0 70000))  
graph export "permanentwage_byage_real.pdf", replace

save "permanentwage_byage.dta", replace
